//
//---------------------  1D  ---------------------
//
//
// Quick is not defined for 1D
//

//
//---------------------  2D  ---------------------
//
//       Staggered Mesh for u-vel and v-vel
//
//       0       1       2       3       4          
//
//5      >       >       >       >       > 
//       |       |       |       |       |              
//   ^---+---^---+---^---+---^---+---^---+---^  4       Mesh for scalar fields
//       |       |       |       |       | 
//4      >   o   >   o   >   o   >   o   >               5  x-x-+-x-+-x-+-x-x
//       |       |       |       |       |               4  x o | o | o | o x
//   ^---+---^---+---^---+---^---+---^---+---^  3           +---+---+---+---+
//       |       |       |       |       |               3  x o | o | o | o x
//3      >   o   >   o   >   o   >   o   >                  +---+---+---+---+
//       |       |       |       |       |               2  x o | o | o | o x
//   ^---+---^---+---^---+---^---+---^---+---^  2           +---+---+---+---+
//       |       |       |       |       |               1  x o | o | o | o x
//2      >   o   >   o   >   o   >   o   >               0  x-x-+-x-+-x-+-x-x
//       |       |       |       |       |                  
//   ^---+---^---+---^---+---^---+---^---+---^  1           0 1   2   3   4 5 
//       |       |       |       |       |                 
//1      >   o   >   o   >   o   >   o   >                 o central node
//       |       |       |       |       |                 x boundary node
//   ^---+---^---+---^---+---^---+---^---+---^  0          > u velocity 
//       |       |       |       |       |                 ^ v velocity
//0      >       >       >       >       >    
//   0       1       2       3       4       5       
//                                                  
//                      v(i,j) 
//                  |     n     |               
//                --|-----^-----|--             
//                  |           |               
//                  |           |               
//     u(i-1,j) = w >     o     > e = u(i,j)
//                  |   (i,j)   |
//                  |           |
//                --|-----^-----|--
//                  |     s     |
//                     v(i,j-1)
//               
//
//---------------------  3D  ---------------------
//
template<class T_number, int Dim>
inline 
bool QuickX<T_number, Dim>::calcCoefficients(const ScalarField &nut) { 
    T_number dyz = dy * dz, dxz = dx * dz, dxy = dx * dy;
    T_number dyz_dx = dyz / dx, dxz_dy = dxz / dy, dxy_dz = dxy / dz;
    T_number ce, cw, cn, cs, cf, cb;
    T_number cem, cep, cwm, cwp, cnm, cnp, csm, csp, cfm, cfp, cbm, cbp;
    T_number nutinter;
    T_number dxyz_dt = dx * dy * dz / dt;

    for (int i =  bi; i <= ei; ++i)
	for (int j = bj; j <= ej; ++j)
	    for (int k = bk; k <= ek; ++k)
	    {
		ce = ( u(i+1, j, k) + u(i,j,k) ) * 0.5 * dyz;
		cw = ( u(i-1, j, k) + u(i,j,k) ) * 0.5 * dyz;
		if ( ce > 0 ) { cem = 0.0; cep = ce * 0.125; }
		else {          cem = -ce * 0.125; cep = 0.0; }
		if ( cw > 0 ) { cwm = 0.0; cwp = cw * 0.125; }
		else {          cwm = -cw * 0.125; cwp = 0.0; }

		cn = ( v(i,j,k) + v(i+1,j,k) ) * 0.5 * dxz;
		cs = ( v(i,j-1,k) + v(i+1,j-1,k) ) * 0.5 * dxz;
		if ( cn > 0 ) { cnm = 0.0; cnp = cn * 0.125; }
		else {          cnm = -cn * 0.125; cnp = 0.0; }		
		if ( cs > 0 ) { csm = 0.0; csp = cs * 0.125; }
		else {          csm = -cs * 0.125; csp = 0.0; }

		cf = ( w(i,j,k) + w(i+1,j,k) ) * 0.5 * dxy;
		cb = ( w(i,j,k-1) + w(i+1,j,k-1) ) * 0.5 * dxy;
		if ( cf > 0 ) { cfm = 0.0; cfp = cf * 0.125; }
		else {          cfm = -cf * 0.125; cfp = 0.0; }		
		if ( cb > 0 ) { cbm = 0.0; cbp = cb * 0.125; }
		else {          cbm = -cb * 0.125; cbp = 0.0; }
//
// nut is calculated on center of volumes, therefore, nut
// must be staggered in x direction:
		nutinter = 0.5 * ( nut(i,j,k) + nut(i+1,j,k) );

		aE (i,j,k) = (Gamma + 2 * nutinter) * dyz_dx 
		    - ce * 0.5 + cep - 2 * cem - cwm;
		aW (i,j,k) = (Gamma + 2 * nutinter) * dyz_dx 
		    + cw * 0.5 + 2 * cwp - cwm + cep;
		aN (i,j,k) = (Gamma + nutinter) * dxz_dy 
		    - cn * 0.5 + cnp - 2 * cnm - csm;
		aS (i,j,k) = (Gamma + nutinter) * dxz_dy 
		    + cs * 0.5 + 2 * csp - csm + cnp;
		aF (i,j,k) = (Gamma + nutinter) * dxy_dz 
		    - cf * 0.5 + cfp - 2 * cfm - cbm;
		aB (i,j,k) = (Gamma + nutinter) * dxy_dz 
		    + cb * 0.5 + 2 * cbp - cbm + cfp;
		aP (i,j,k) = aE (i,j,k) + aW (i,j,k) + 
		             aN (i,j,k) + aS (i,j,k) +
		             aF (i,j,k) + aB (i,j,k) + dxyz_dt 
		             + cem - cwp + cnm - csp + cfm - cbp;
//		aP (i,j,k) /= alpha;  // under-relaxation
//		+ (ce - cw) + (cn - cs) + (cf - cb);
		sp (i,j,k) = u(i,j,k) * dxyz_dt -
		    ( p(i+1,j,k) - p(i,j,k) ) * dyz +
		    nutinter * ( (v(i+1,j,k) - v(i,j,k) - 
				  v(i+1,j-1,k) + v(i,j-1,k)) * dz +
				 (w(i+1,j,k) - w(i,j,k) - 
				  w(i+1,j,k-1) + w(i,j,k-1)) * dy ) +
		    u(i,j,k) * (1-alpha) * aP(i,j,k)/alpha;// under-relaxation

		if (i < ei-1)        { sp (i,j,k) += cem * u(i+2,j,k); }
		else if ( i == ei-1) { sp (i,j,k) += cem * u(i+1,j,k); }
		if (i > bi+1)        { sp (i,j,k) -= cwp * u(i-2,j,k); }
		else if ( i == bi+1) { sp (i,j,k) -= cwp * u(i-1,j,k); }

		if (j < ej-1)        { sp (i,j,k) += cnm * u(i,j+2,k); }
		else if ( j == ej-1) { sp (i,j,k) += cnm * u(i,j+1,k); }
		if (j > bj+1)        { sp (i,j,k) -= csp * u(i,j-2,k); }
		else if ( j == bj+1) { sp (i,j,k) -= csp * u(i,j-1,k); }

		if (k < ek-1)        { sp (i,j,k) += cfm * u(i,j,k+2); }
		else if ( k == ek-1) { sp (i,j,k) += cfm * u(i,j,k+1); }
		if (k > bk+1)        { sp (i,j,k) -= cbp * u(i,j,k-2); }
		else if ( k == bk+1) { sp (i,j,k) -= cbp * u(i,j,k-1); }
	    }

    calc_du_3D();
    applyBoundaryConditions();

    
    return 1;
}
















